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SUMMARY 

An approximate factorization scheme based on the AF2 algorithm is presented for solving the three- 
dimensional full potential equation for the transonic flow about isolated wings. Two spatial 
discretization variations are presented; one using a hybrid first-order/second-order-accurate scheme 
and the second using a fully second-order-accurate scheme. The present algorithm utilizes a C-H grid 
topology to map the flow field about the wing. One version of the AF2 iteration scheme is used on 
the upper wing surface and another slightly modified version is used on the lower surface. These two 
algorithm variations are then connected at the wing leading edge using a local iteration technique. 
The resulting scheme has improved linear stability characteristics and improved time-like damping 
characteristics relative to previous implementations of the AF2 algorithm. The presentation is 
highlighted with a grid refinement study and a number of numerical results. 


INTRODUCTION 

The lonq term objective of this research is to develop a chimera-based full potential flow solver which 
will be compatible with the well-established OVERFLOW Euler/Navier-Stokes flow solver developed 
at NASA Ames. 1 Thus, the user will have an option of which flow solver to use in the chimera-based 
zonal grid approach: full potential, Euler or Navier-Stokes. Of course, the full potential option will not 
be applicable for all applications, but for those that are applicable, the execution time should be up to 
two orders of magnitude less than for the Navier-Stokes formulation. Indeed, a chimera-based full 
potential solver should have modest execution times on even moderate-speed workstations. In a 
parametric study the bulk of the required computations could utilize the full potential approach and 
then a few selected conditions could be "checked' 1 with a more complete, and thus more accurate, 
Euler or Navier-Stokes simulation. Such an approach would be extremely cost effective especially 
considering that all of these approaches would utilize the same problem setup and post processing 
software and to a large extent the same grid generation software. Applications of this new approach 
are quite numerous and include providing a fast mechanism for assessing wind tunnel wall and 
support interference effects associated with wind tunnel testing, or it could be used directly in the 
industrial preliminary design environment. 


The specific goal of this report is to document recent advances in approximate factorization algorithms 
for solving the full potential equation that will eventually be useful in a multi-zone chimera 
environment. However, only single-grid-zone methodology will be presented herein. In particular, 
the presentation will focus on numerical solution of the transonic flow about an isolated wing utilizing a 
sinqle-zone C-H-topology grid. For this grid all constant span stations on the wing exhibit the familiar 
C-qrid topology, and all constant chord stations across the wing exhibit the equally familiar H-grid 
topology. Utilization of the chordwise C-grid topology is useful because it lends itself more readily to 
boundary layer correction implementation, which is future goal of this work. 

In the present study the AF2 full potential algorithm developed in Refs. 2-3, which uses an O-H grid 
topology for isolated wing computations, is modified for the present C-H grid topology. Direct 
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application of the Ref. 2-3 AF2 iteration scheme on C-type grids is difficult because of a cell-aspect- 
ratio stability limitation, which exists at the wing surface and is greatly accentuated by the C-type grid 
topology. The nature of this instability, which was first described by South and Hafez, 4 will be briefly 
described in the iteration scheme section of this report. The present AF2 scheme variation is 
designed to control this instability and provides efficient and reliable convergence for a wide range of 
isolated-wing transonic flow simulations involving C-H topology grids. 

This presentation begins with a discussion of the governing equation formulation followed by a 
detailed discussion of the numerical algorithm including spatial discretization scheme, boundary 
conditions, vortex sheet conditions, and the newly modified AF2 iteration scheme. Next, typical 
transonic wing computational results simulating the flow about the ONERA M6 Wing are presented. 
These results include a grid refinement study showing the levels of error in lift and drag relative to 
available Euler results. Finally, the presentation ends with concluding remarks and recommendations 
for future work. 


GOVERNING EQUATION FORMULATION 

The steady, three-dimensional, full potential equation written in strong conservation-law form is given 
by 


{P<P X ) X +(P<Py)y+ (P<I> Z ) Z =0 


P = 


1 - 


+ <Py + <Pz) 


1 

r-i 


(la) 


(1b) 


where p is the fluid density; x, y, and zare Cartesian coordinates; y is the ratio of specific heats; and 

— ^ 

<t> is the full or exact velocity potential related to the velocity vector q by 

V(p = q 

The velocity components can be expressed using Cartesian coordinates as follows: 


V<p = q = ui + vj + wk = <t> x i + <t> y j+<p z k 


where /, / .and k are the standard unit vectors in the x, y, and z directions, respectively. The mere 
existence of the velocity potential implies that the curl of the velocity vector must vanish. Thus, flows 
governed by the full potential equation must be irrotational. In addition, derivation of the density 
relation [Eq. (1b) above] requires the assumption of isentropic flow. In Eqs. (1) the density (p) and 
velocity components ( <p x , <p y , and <p z ) are nondimensionalized by the stagnation density (p s ) and the 
critical speed of sound (a.), respectively. Additional relations valid for these flow assumptions and 
this nondimensionalization include 

Isentropic equation of state: 


p r+i 

P r 2 y 
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Bernoulli's equation: 


q 2 ! a 2 1 y+1 
2 + y-1 _ 2 y-1 


Speed of sound definition: 

2 P 
a = y— 

P 

where p is the fluid pressure, a is the local speed of sound, and q is the magnitude of the local fluid 
velocity. 

To complete the above system, boundary conditions are required along all boundaries. Specifically, 
these boundaries fall into three categories: freestream, symmetry planes, and geometric surfaces. 
The freestream boundary condition, simply stated, is given by 

x 2 + y 2 + z 2 0 -» K 


where <p„ is the freestream distribution of the velocity potential, usually uniform flow. The latter two 
boundary conditions, symmetry planes and geometric surfaces, are both treated in the same manner, 
i.e., with a flow tangency assumption given by 


q* n = 0 


where n is a unit vector normal to the geometry of interest. More on boundary conditions including 
numerical implementation will be presented in the Numerical Approach section. 

Equations (1) express mass conservation for flows that are isentropic and irrotational. Despite these 
limiting assumptions, the full potential formulation can be used in a shock-capturing context providing 
the shock waves are weak. The corresponding shock-jump conditions are valid approximations to the 
Rankine-Hugoniot shock jump conditions (derived from the Euler equations) for many applications. 
The key parameter in this situation is the normal component of the Mach number just upstream of the 
shock wave in question, which must remain below about 1.3 for the full potential formulation to be a 
reasonable approximation to the Euler equations. This is well within the scope of many transonic flow 
applications and includes the cruise conditions for most transonic transport aircraft. More discussion 
on this point including a comparison of the Euler and isentropic full potential shock polars is 
presented in Steger and Baldwin. 5 

Equations (1) are transformed from the physical domain (Cartesian coordinates) into the 
computational domain using a general, independent-variable transformation. This general 
transformation, indicated by 


$ = $(x,y,z) 

q = q(x,y,z) (2) 

£ = C(x,y,z) 

maintains the strong conservation-law form of Eqs. (1). The final transformed version of the full 
potential equation is given by 


pU 

J 



= 0 


(3a) 
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where the density expression becomes 

i 


P = 




and 


U = + /t 5 0£ 

V = A4<p£ + A 2 <t>ri + Ae^ 

W = a 5<P% + A 6 < t > j] + a 3 

= + ^ + « | 

*2 = V?j»Vtj=7 ll + Vy + ilz 

43»Vf.Vf-g + ?J + f| 
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~^z r lyCx ~^y T lxiz ~ l ^x r lzCy 
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(3b) 
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In Eqs. (4), U, V, and W are the contravariant velocity components along the <£, 77, and f coordinate 
directions, respectively; A, - Ag are metric quantities; and J is the determinant of the transformation 
Jacobian. The above full potential governing equation formulation can be used for general 
geometries in which the aerodynamic surface of interest is mapped to a constant coordinate line in the 
computational domain. This makes the flow-tangency boundary condition easy and accurate to 
implement. In transformed coordinates the flow-tangency boundary condition becomes (e.g., for an 
rj=constant surface) 


q» V77 = 0 

-> ( -* -» -O 

q»Vj 1 = U x i + <f> y j + <p z k 

V = 0 



i + r] y i+n z k 


= 1/ 


More simply stated, the contravariant velocity component in the rj -direction V must vanish at the 
^constant surface where flow tangency is required. The above flow tangency boundary condition is 
also used at y= 0 as a symmetry plane boundary condition. 

For all computations presented in this report the % coordinate is aligned with the C-grid wrap-around 
direction. Hence, it is aligned with the negative flow direction below the wing and with the positive 
flow direction above the wing. The 77 coordinate is aligned with the spanwise direction, and the £ 
coordinate is in the wing normal-like direction. The orientation of each of these coordinate directions, 
especially the £ coordinate, will have an influence on the construction of the iteration scheme as will 
be seen in the next section. 
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NUMERICAL APPROACH 


Spatial differencing scheme 

A second-order-accurate spatial differencing approximation to the full potential equation written in the 
general transformed coordinate system is given by 


■m 


+ 8, 


+V2,j,k 


pv) 

J )i,j 


+ S( 


7+1/2, V 


'pW_\ 

. j )i,j, 


= 0 


(5) 


,k+V2 


where (for example) 

Ui+V2j.k = A ( < Pi+\j,k “ Qi.j.k) 

+ ~ A, +VJ/> (Qi+Xi+Xk ~ Pi+Xi-Xk + Qi.j+Xk ~ Qi.j-Xk) 
+~ A /+V 2.m ^'+1./'.* +1 ~ Qi+Xi.k-1 + ^i.j.k+l ~ Qi.j.k-l) 


(6) 


The operators are standard backward difference operators in the three coordinate 

directions defined by 


<M-),7. 

*=(—), 7, 

k ~(—)i-\j,k 


S n Uij, 

*=(-)/,/, 

k ~ 

(7) 

S(Uu. 

*=(-)/,/, 

k 1 



The /,/, and k subscripts used in the Eqs. (5)-(7) are used to denote position in the finite-difference 
grid such that £ = /A£, rj = /A 77, £ = kAC,. For convenience, the , At? and A£ values are taken to be 

unity. 

Values of p, and J are computed using the freestream preserving spatial differencing 

scheme described by Flores et al ® and Thomas and Holst. 7 This scheme produces a zero residual for 
each interior grid cell with a freestream distribution of the velocity potential and generally produces a 
solution with improved accuracy, especially near grid singularities or in regions of rapid grid stretching. 

At supersonic points the centrally-differenced, second-order-accurate spatial discretization scheme 
just presented must be upwind biased to maintain stable operation. This is accomplished by replacing 
the density in the £-flux computation with an upwind biased value of the density. Thus, the new 

spatial differencing scheme is given by 






+ 8 , 


'i+V2,j,k 


pV] 

J 


Jij 


+ 8t; 


7+1/2,/c 


pW) 

J 


hi. 


= 0 


i,k+V2 


(8) 


where the density coefficient p, +V2 . M is defined by one of two options. The first option is given by 


Pi+V2,j,k -Pi+V2,j,k v i+V2,j,k(Pi+V2,i,k Pi-V2,j,k ) 


(9) 


5 



where 


_ |2.46625(2p P/+ 1 / 2 ,/,* Pi-V 2 ,j,k)^ ^ Pi,i,k—P Mfh 

Vw ' 2 '« 1 0 if K,y>9 (0 ’ 

The quantity p* is a constant equal to the sonic value of the density, which for y = 1.4, is 
0.6339382 -.. The quantity C is a user specified coefficient usually set near a value of one, and 
V i+V 2 jk is a switching parameter that controls the amount of upwinding that exists in the numerical 

scheme. Equation (10) is designed such that the value of v /+1/2 ,/,* will be zero at all subsonic grid 

points, i.e., a second-order-accurate, centrally-differenced scheme is retained for subsonic regions of 
the flow domain, and larger than zero at all supersonic grid points, i.e., a first-order-accurate upwind 
scheme is utilized for supersonic regions of flow. The supersonic branch of Eq. (10) approximates 
-1 )C, and thus, the amount of upwinding increases dramatically as the extent of supersonic 

flow increases. To keep the value of p /+1 / 2 ./> bounded by p i+V 2,i,k and the value of ^+1/2, y.* 

is constrained to be less than or equal to one. 

The second option for defining the density coefficient (p, + 1 / 2 , 7 ,*) (inspired by Kinney et al . 8 and 
Jameson 9 ), is given by 


P/+1/2,/,* ~ Pi+V2,i,k ~ v i+V2,j,k[Pi+V2,i,k Pi-M2,j.k 

~ ^ /+1/2, /, * (Pi-1/2,/,* ~ Pi-2/2,j,k )1 


where the limiter 4* is defined by 


4 1 


(+ 1 / 2 ,/,* 


jl - C 2 A if r 7+1/2, 7, k - 0 

| 0 if r ,-+1/2,/,* ® 


and r i+V2 j tk is the ratio of successive density gradients defined by 


r (+1/2,7,* 


P/+1/2,/,* ~ P/-1/2,/,* 

P/-1/2,/,* ~ P/-3/2.7,* 


( 11 ) 


(12) 


(13) 


In Eq. (1 2) the quantity C 2 is a constant typically set to a value near one and A is the local £ -direction 
grid spacing which is approximated using 


A = (ITEU - ILE)-' 

where ILE is the /th grid index at the wing leading edge and ITEU is the /th grid index at the wing 
upper-surface trailing edge. To improve stability, especially for finer grids, the limiter function is 
decreased in magnitude with increasing distance away from the wing surface. A function of the form 

y i+V2,j,k = ^(+1/2, 7, 1 (^ 3 ) <k 11 04) 

where C 3 is a constant set to a value just above one, e.g., 1.08, and k is the normal-like £ -direction 
grid index with k= 1 defining the wing surface, has worked well for this purpose. The resulting scheme 
retains second-order-accuracy at the wing surface while allowing increased stability associated with 
increased dissipation away from the wing surface. 
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At subsonic grid points, the second density coefficient option (identical to the first option) produces a 
zero value of v i+V2J k , leading to a second-order-accurate, centrally-differenced, spatial-differencing 

scheme. At most supersonic points rwill be greater than zero and the resulting spatial-difference 
scheme will be upwind-biased and second-order-accurate. At a supersonic point which is an 
extremem r<0. This produces a zero limiter function {*¥ = 0) at these points, which in turn produces 
the original first-order density upwinding option. 

The lead truncation error term generated by the first density upwinding option at supersonic grid 
points can be approximately written as 




This expression is characteristically dissipative and leads to a first-order-accurate scheme. Larger 
values of this term, achieved by increasing the value of C in Eq. (10), will produce increased smearing 
of the solution in supersonic regions of flow. 

The lead truncation error term generated by the second density upwinding option at supersonic grid 
points can be approximately written as 

A£ 2 (vp^)£ + A£C 2 A(vp^)£ 

The first term is characteristically dispersive, and the second term is characteristically dissipative. 
However, providing C 2 A approaches zero as approaches zero, both terms are second-order 

terms and the resulting scheme is second-order accurate. 

The two supersonic spatial differencing schemes represented by Eqs. (8)-(10) and Eqs. (8), (10)-(14) 
are valid for supersonic flows which are approximately aligned with the positive £ -coordinate 
direction. In practice, even for C-H topology wing grids involving large amounts of sweep, this type of 
supersonic upwinding is suitable. It is the only type of supersonic flow stabilization used in the 
present study. Nevertheless, generalization of the present scheme for arbitrary orientations of a 
curvilinear coordinate system is easy to accomplish. See Refs. 2-3 for examples involving a 
generalized form of the first-order density upwinding option. 


Boundary conditions 

Flow tangency and symmetry plane boundary conditions (as described above) require that the 
velocity component normal to the applicable boundary must vanish. For a general nonorthogonal 
mapping, such as that described by Eq. (2), the general condition for flow tangency requires the £ 
contravariant velocity component at C = Cmin (for example) be zero (i.e., =0). This is 

implemented in the present study using a mass-flux reflection condition given by 




(15) 


where in this case /c=1 corresponds to the tangency or symmetry plane surface. In other expressions 
where <p ( is required at /c=1 [e.g., in Eq. (6) or in the density computation at k= 1], the Wr =(nm =0 

boundary condition is used again to obtain 



( 16 ) 
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Thus, a value of <p ( at the boundary condition surface can be obtained without using one-sided 
differences on the velocity potential. 

At the wing-surface/symmetry-plane line of intersection, two contravariant velocity components must 
vanish (i.e., V v=v ^ =0, W f=fmin =0). This is implemented using two mass-flux reflection conditions 

given by 


pV 

J 






(17) 


where in this case /= 1, k= 1 corresponds to the line of intersection. In other expressions where 
and are required at the line of intersection [e.g., in Eq. (6) or in the density computation at j=1, 
k=1 ], the V v=Vwin = 0, = 0 boundary conditions are used again to obtain 


<P 


77 (intersection 


a z a 3 - a| ^ 


*<\ 


intersection 


~ AgA^ 

a 2 a 3 - a| s 


(18) 


Thus, values of <j> n and at the symmetry-plane/wing-surface line of intersection can be obtained 
without using one-sided differences on the velocity potential. 

At the outer boundary a standard freestream boundary condition is used which is given by 

<t>(x,y,z) = u^x + w„z 


where u„ and iv„ are standard Cartesian velocity components in the freestream associated with the 
x and z directions, respectively. 


Vortex sheet conditions 

For lifting computations involving the full potential formulation, circulation is accommodated with the 
usual vortex sheet or wake cut emanating downstream of the lifting surface (for wings, downstream of 
the wing trailing edge). The amount of lift or circulation is equal to the jump in velocity potential across 
the vortex sheet. The jump or discontinuity in velocity potential, as well as the double-stored 
characteristic of the vortex sheet, must be accounted for in the residual mass flux and density 
computational logic described above. For example, the U contravariant velocity component 
computation of Eq. (6) must be modified as follows at the vortex sheet 

^/+ 1 / 2,/,1 = A „ V2;1 ( 0 /+ 1./.1 - 

+ J ^4,. vw (0;+1,/+1,1 ~ Qi+Xj-V + 0/J+1.1 ~ Qi.j- 1.1 ) (1 9 ) 

+ 4 ( 0/+U2 - 0AH-/J.2 - Ty + 2 - <t>MI-i+\j'2 “ Ty ) 
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where the k = k min (k = 1) boundary has been mapped to the wing/vortex-sheet surface using a C- 
qrid topology. The point i,j,1 is atypical point on the upper side of the vortex sheet and Nl-i+1,j,1 is 
the identical point on the lower side of the vortex sheet. The quantity r y is the /th value of the 

circulation (where / is measured along the span of the wing) and i=NI corresponds to the maximum 
value of the /th coordinate. The above logic requires identical grid point distributions to exist on both 
the upper and lower wing/vortex-sheet surfaces with a grid point at the wing leading edge, thus, the 
/th point along the upper (or lower) vortex sheet exactly matches the (Nl-i+1)\h point along the lower 
(or upper) vortex sheet. The value of r y is computed using 

r y = QrrEu.jA - Qitelp 

where ITEU,j,1 corresponds to the /th location of the upper wing trailing edge and ITELJ.1 
corresponds to the same point on the lower wing trailing edge, i.e., ITEU=NI-ITEL+1. The /th-location 
jump in velocity potential established at the wing trailing edge is maintained downstream along the 
vortex sheet, i.e., 


i=iteu+vteu+2,....ni Pi) 

Because of Eq. (21), only one velocity potential residual computation [using Eq. (8)] is required for 
each two grid nodes on the vortex sheet downstream of the wing trailing edge. This single residual 
value for each two vortex sheet points is computed by performing normal residual computations for 
both points and then averaging the values. 


Iteration scheme 

The iteration scheme utilized in the present study is called the AF2 scheme and was first introduced 
by Ballhaus and Steger 10 in 1975 for solving the transonic small-disturbance potential equation for 
two-dimensional applications^ The present implementation is closely related to the AF2 scheme 
described in Holst and Thomas, 3 which was designed for solving the three-dimensional full potential 
equation for isolated-wing applications using an O-H grid topology (as mentioned in the Introduction). 
The present algorithm is also designed for solving isolated wing applications, but uses a C-H grid 
topology. The present implementation’s C-type grid topology is more amenable to boundary layer 
correction implementation and is the major motivation for the present modification. More on the 
differences between the Ref. 2-3 AF2 scheme and the present version will be presented when the 
actual scheme is discussed. 

A general iteration scheme for solving the full potential equation can be expressed as 

NC?j k + (oL$j k = 0 ( 22 ) 

where the n superscript is an iteration index, L4? JM is the nth iteration residual at the position in 

the finite-difference mesh [spatial difference scheme for the full potential equation defined by Eq. 
(8)], L is the residual operator, to is a relaxation parameter, C”j k is the correction defined by 


§ The name AF2 (short for Approximate Factorization Scheme 2) was given to this scheme in Ref. 10 
because it was the second scheme presented in that study. 
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and N is the "left-hand-side" operator that determines the type of iteration scheme. The AF2 iteration 
scheme used in the present formulation (one of two variations) can be specified using the following 
definition for N: 


a 2 NC° jM = - a-S 6 R, 


cx — S rf Sj S 7 j 


a-Sa k 8^\-a 2 El' 


CPjjc (23) 


where a is an acceleration parameter (to be discussed shortly); R,, Sj, and T k are coefficients 
defined by 





S l = 



i,j-V2,k 


T k = 



n 


(24) 


5$ , Stj, are standard forward-difference operators in the three coordinate directions defined by 


SiUij, 

k -(-)i+\j,k ~ (-)/ 

8nU U . 

k ~ (— )i,j+\k ~ (— )i 


k = (_)/,/,* +i “(_)/ 


and is a shift operator defined by 


(25) 


A standard linearized von Neumann stability analysis of the interior AF2 scheme given by Eqs. (8), 
(22)-(23) yields a stable scheme providing a>0 and 0<tu<2. The quantity a behaves like the 
inverse of a time step in a typical time-dependent iteration, this implies that the linearized AF2 scheme 
is stable for any positive time step and hence is said to have unconditional linear stability. 

The fact that a behaves like the inverse of the time step means that small values of a (large time 
steps) advance the solution rapidly and often cause (nonlinear) high frequency error growth. Large 
values of a (small time steps) advance the solution slowly and provide solution smoothing especially 
for the high frequency error components. Thus, a sequence of values for a produces optimal 
steady state convergence and can be obtained using 


«* = «H 


a. 


k - 1 
M - 1 


V a H J 


k = \2,--,M 


(26) 


where a L and a H are low- and high-frequency limits for the a k parameter sequence, respectively, 
and M is the number of elements used in the a sequence (see Refs. 2-3 and 1 1 for more details and 
a variety of applications). 

The AF2 factorization given by Eq. (23) is implemented in a three-sweep format, each involving a set 
of banded, scalar matrix inversions. These sweeps are given by 
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Sweep 1: 


(27a) 




Sweep 2: 

(a — S t] Sj S n )g n jM — h,j,k + a 

(27b) 

Sweep 3: 

(a- 8{T k 8 f )Cf jk =g n jM 

(27c) 


In step 1, the f array is obtained from the residual by solving a simple bidiagonal matrix equation for 
each £, line. The g array is then obtained from the f array by solving a tridiagonal matrix equation for 
each tj line. Finally, the correction array is obtained in the third sweep from the g array by solving a 
tridiagonal matrix equation for each £ line. Note that sweeps 2 and 3 can be performed together, i.e., 
immediately after the g values for the Ah computational plane are obtained from sweep 2, values for 
the correction array can be obtained for the same Ah computational plane from sweep 3. Thus, the g 
intermediate result array needs to be only two dimensional. 

As can be seen from Eq. (23) or Eqs. (27) the left-hand-side difference operator in the £ direction 
(the C-grid wrap-around direction in the present implementation) has been split between the first and 
second factors. This operator splitting is a telltale characteristic of the AF2 scheme. However, the 
operator splitting does not have to be in the first factor; any of the factors can be split. For example, 
the AF2 scheme of Ref. 3 splits the £ direction factor. Like the present approach, the Ref. 3 
implementation uses the £ direction as the wrap-around direction, rj is spanwise, and £ is in the wing 
normal-like direction. However, unlike the present implementation, the Ref. 3 approach used an O-H 
grid topology. 

Implementation of the AF2 scheme (like all implicit schemes) requires matrix "boundary conditions" for 
each sweep. In particular, conditions are required for f at / = A//, g at j = 1 , g at j = NJ ; C at / = 1 , C at k = 1 
and Cat k= NK. The second and third of these conditions are associated with sweep 2 and are 
satisfied using a standard Dirchlet condition if the boundary is freestream or a standard Neumann 
condition if the boundary involves flow symmetry or flow tangency. The fifth and sixth conditions are 
associated with sweep 3 and likewise, are satisfied using a Dirchlet condition if the boundary is 
freestream or a Neumann condition if the boundary involves flow symmetry or flow tangency. The first 
and fourth conditions are nonstandard and result from the special AF2 left-hand-side operator 
splitting. These two special left-hand-side matrix "boundary conditions" are discussed next. 

First, in sweep 1 , a value of fat l=NI is required for each value of j and k. A simple technique to satisfy 
this condition is given by 


SiR,f 


n 


= 0 


This condition works, but a stability limitation at the f=N/ boundary results, which is of the form 4 


a>Rji o 

where R ) is defined by Eq. (24) and (o is the relaxation factor from Eq. (22). Thus, there is a limit on 
the size of a, which can have a dramatic effect on convergence (even if the limitation exists only at 
the computational boundary). The key aspect of this stability limit is the A, / J contribution inside R,, 
which behaves like a grid-cell aspect ratio (boundary tangential spacing over boundary normal 
spacing). 
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In the Ref. 3 AF2 implementation with the wing normal-like f operator split, the boundary stability 
limitation described above, exists at the wing surface. For the O-type grid in use for the Ref. 3 
application, in which the cell aspect ratio is 0(1) all around the inner boundary, this stability limitation 
has a minor consequence. However, direct application of the Ref. 3 AF2 iteration scheme for solving 
the full potential equation on C-type grids would be very inefficient because the cell aspect ratio for 
the inner boundary of a C-grid becomes quite large downstream of the wing trailing edge along the 
wake cut. 


In the present approach the AF2 scheme is configured such that this approximate matrix inversion 
boundary condition is implemented away from the wing surface. In particular, it is implemented along 
the downstream outflow boundary where the flow is freestream and the cell aspect ratio is generally 
quite small. 

The second special "boundary condition" for the left-hand side is associated with sweep 2. A value 
for the correction at /= 1 is required. If / = 1 represents an outer freestream boundary, then this 
condition is succinctly satisfied without loss of stability by using C= 0, i.e., the solution is required to 
be freestream along the /' = 1 boundary with corrections for all iterations being zero. For the present 
AF2 implementation, as will be seen shortly, this condition must be satisfied along a surface of grid 
points that emanates from the wing leading edge that exactly divides the C-grid topology into two 
halves. The intermediate condition on the correction at this location used in the present scheme, 
involves a local iteration to ensure implicitness and will be subsequently discussed in detail. 

An additional special condition is required at the vortex sheet downstream of the wing trailing edge. 
Values of the correction array are discontinuous across the vortex sheet during the iteration process 
due to the changing value of the circulation, namely 


C n 

ij, 1 ” 


pH + 1 T'lff 

- 1 / 1 / 


where k= 1 corresponds to the wing/vortex-sheet-cut grid surface and T" is the nth iterate of the 

circulation at the yth spanwise station. In the above expression the i subscript must be taken along the 
upper vortex sheet, and thus, A//-H-1 automatically is along the lower vortex sheet. This discontinuity 
must be taken into account when the sweep 3 matrix coefficients are computed downstream of the 
wing trailing edge. 

Some comments about the AF2 scheme given by Eq. (23) are in order. First of all, the role played by 
each of the terms in this factorization can be understood more clearly by multiplying out each of the 
factors, which yields 

a 2 NC? ]k = a 3 - a 2 8 n S t 5 - a 2 8 ( T k 8 ( + a 8 „ S J 8 n S ( T k 8 \ 

-a 3 E^ - a 2 5$Rj + a8sRi5 n Sj8 n +a8‘Ri8sT k 8i; (28) 

- 8 s R, 8 „ S, 8 „ S { T k 8 ( + a 2 8$ R,Ej' )c£- k 


Terms 1 and 5 when combined become 

a 3 - a 2 Ef = a 3 (l - E^) = a 3 8$ 
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which when operating on the correction produces a 0*,- like term that provides implicit (i.e., left-hand- 

side) or time-like dissipation to the iteration process in supersonic regions of flow § Existence of this 
term in the AF2 factorization is the chief reason that the AF2 scheme enjoys a convergence efficiency 
advantage over more traditional ADI-type algorithms for solving the full potential equation for transonic 
flows with shock waves. This term provides favorable implicit dissipation in supersonic regions of 
flow providing the S, coordinate is approximately aligned with the local streamwise direction “ s " and 
providing the backward difference corresponds to an upwind difference. In other words, the <p it - term 

must approximate an upwind-evaluated 0 s( -term in supersonic regions of the flow. See Jameson 13 
for more information on this point. 

The 4> it term arising from the Eq. (23) factorization is a direct result of having the £ -direction factor 
split between sweeps 1 and 2. In contrast, the Ref. 3 AF2 scheme had the £ -direction difference 
split between these two sweeps. Thus, the <j>^ t term (or equivalently the <p st - term) required for a 

convergent iteration in the Ref. 3 AF2 scheme had to be added as an extra term to the appropriate 
left-hand-side factor when the flow became supersonic. 

In summary, the present scheme's £ direction operator splitting is more efficient for two reasons. 
First, it requires fewer operations because no additional time-like dissipation teims have to be added 
in the if; -direction. Second and most important, the boundary stability limitation resulting from the 
AF2 operator splitting is effectively removed. The present AF2 factorization with a <!; -direction 
splitting, represents a significant improvement over past AF2 approaches for solving the full potential 
equation on C-topology grids. 

Terms 6 and 10 from Eq. (28) become 


-a 2 8$ Rj + a 2 8$ R t E J 1 = -a 2 8 $ Rj{A-E^ = -a 2 8$ Ft, 8$ 


which, assuming the cross derivative metrics A 4 and A 5 are small, is a close approximation to the first- 
difference-operator term of Eq. (8). In the limit as an orthogonal grid is approached, i.e., as A 4 and A§ 
go to zero, the approximation becomes exact. 

This term coupled with terms 2 and 3 from Eq. (28) provides an approximation to the entire right-hand- 
side residual operator and can be written as 


- a 2 8 $ R, <5 - + 8 Sj 8 n + 8 c;T k 8^ C” / k ~ -a 2 LC t j k 


where a strict equality exists in the special case of an orthogonal grid. 

The remaining terms in Eq. (28) (terms 4, 7, 8, and 9) are factorization error terms (FET) which are 
driven to zero as the iteration proceeds. Thus, Eq. (23) can be written as 


a 2 NC^ k - a 2 LC°j' k -a 3 8 ^ + FET 


(29) 


§ The time-derivative in the <p,, term is obtained by assuming the relaxation or iteration process to be an 
iteration in real time. Thus, 


a 8$ cPjjt = A t - 4>u.k + 4>i'-\i.k ) ~ Qy 
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which shows that the AF2 scheme left-hand-side N operator represents a close approximation to the 
right-hand-side residual operator. This is an indication of rapid convergence. 

Because of the inherent direction on the <p it term differencing operator, i.e., always backwards, the 
AF2 scheme given by Eq. (23) is suitable in supersonic regions of flow only when the positive £ 
direction and the flow direction are approximately aligned. For supersonic regions of flow where the 
negative % direction is approximately aligned with the flow direction, i.e., below the wing, the AF2 
scheme given by 


a NC U,k = - a+6{R M 


a-SnSj T k 5fj-a 


>E? 


C?jjc (30) 


is an appropriate alternative. Note that the only changes between Eq. (30) above and Eq. (23) are in 
the first factor and in the sign of the shift operator superscript. When the factors of Eq. (30) are 
multiplied out the resulting expression for the N operator is given by 


a 2 NCf j k = a 2 LCf j k + a 3 5$ Cf j k + FET (31 ) 

Thus, the same level of approximation to the residual operator is obtained by the factorization of Eq. 
(30). The only significant difference is in the direction and sign of the <p^ t term which is now upwind 

(as desired) for all regions of flow in which the negative £ coordinate is aligned with the positive flow 
direction. 

The AF2 factorization of Eq. (23) is appropriate (for example) for solving the transonic flow on the 
upper surface of a wing using a "C" or "O" type grid topology, and the factorization given by Eq. (30) is 
appropriate for the lower surface. This assumes the grid lines are wrapped in a clockwise direction 
with the flow from left to right. The two schemes need only be "connected” to each other to allow for 
general wing (or for that matter, any lifting surface) computations. This is accomplished using the 
following algorithm written in a three-step format, which is intended for C-type grid topology 
applications about isolated wing geometries. The / subscript is assumed to be in the wrap-around 
direction, j is assumed to be spanwise, and k is assumed to be in the normal-like direction. Each of the 
indicated operations is to be performed for all values of j and k, but only for the values of / that are 
indicated. 


Step 1 (Sweep 1): 


(a-S.R,)^ 

M = a2a)L< Pi,i.k 

/ = A// — 1.A//-2, -,/LE 

(cc+Si; R M )f' u 

,k = a 2 a)L<l>” jk 

/ = 2, 3, -.ILE 


f n _ fU 
T i,j,k “ r i,j,k 

/ = A//-1.A//-2, -,/LE-l 



/ = 2,3, -,/LE -1 

Kjm = 


/ = ILE 
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Step 2 (Sweep 2/3 leading edge iteration, m=1,2,...,MAXIT l ): 


(a-S n S ) 8 n )g“ ; r=f? J>k + a 2 CWj i l 
(a- S v Sj 5 n )g‘j™ = f" jk + a 2 C?;™~l 
(a-8sT k S ( )C“f k =g1:r 
(a-8 i T k 8 ( )C l ^ k =g^ 

c,T*=5( c T« +c S) 


i = ILE 


C n m 

i,j,k ~ ° ij,k 

C*=</.* +c "m 


/ = ILE, if m = MAXIT 


(a -S n Sj 8 n )g u h ™ = f"j,k + « 2 % 


(a-S<T k S ( )C?f k = g 


i.K 


i = ILE + 1 


(a- 5,S y 5 n )flr^=f / '; M + a 2 C^ 


7+i.y.ir 

(a-8 ( T k 8 ( )C?£ k =gtf 


/ = ILE - 1 


(32b) 


Cn _ s^n,m 

i,j,k ~ ° i,j,k 

<t>ul = <PiJ.x +C U.x 


i = ILE + 1, ILE - 1, if m = MAXIT 


Step 3 (Sweep 2/3 downstream of the leading edge): 

(a - 8 „ S y 5 „ )g u jk = + a 2 C," t y * 

(a-8;T k 8 ( )C?j ik =gl k • 

= 0/J./r +C U,k 
{a-8 n Sj8 n )g l jk = f"y >k + a 2 C/+i,y,* 


i = ILE + 2, ILE + 3,-,/V/ -1 


(a- 8t;T k 8z)Cj , Jk = g‘ jjk 



+ C”j, 


/.* 


i = ILE -2, ILE -3, -,2 


(32c) 


where i=ILE corresponds to the wing leading edge, /=M corresponds to the maximum value of /, the 
superscripts uand / denote upper- and lower-wing-surface values, and the superscript mis an index 
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for the local iteration used to obtain the nth level solution at i = ILE,ILE± 1 . The initial values of the 
correction at the leading edge used in the local iteration are simply assumed to be zero, i.e., 


The quantity MAXIT is the maximum number of iterations used for the local leading edge iteration. A 
value between 10 and 20 is typically used. Use of larger values will not improve the global 
convergence rate. 

In step 1 [Eq. (24a)], the /array is obtained by solving a simple bidiagonal matrix equation for each £ 
line. The upper and lower surface computations are completely independent of each other. The two 
values of / obtained along the grid surface emanating from the wing leading edge (i=ILE,j,k) are 
averaged to produce a single unique value of f at this location. Next, as summarized in step 2, the 
correction array is obtained by local iteration at the wing leading edge.§ This is achieved by alternate 
implementations of sweeps 2 and 3 at / = ILE, ILE ± 1. After the leading edge correction is obtained, 
the remaining corrections in the three-dimensional field are obtained by sweeping the sweep 2/3 
combination away from the leading edge on both the upper and lower wing surfaces. At the wing 
trailing edge the upper and lower sweep 3 matrix inversions are combined into a single inversion, thus 
providing the maximum amount of implicitness for the iteration scheme downstream of the trailing 
edge. 

Despite the improved time-like dissipation arrangement of the scheme just presented, additional 
dissipation is required for some cases. For swept wing cases the £ -coordinate direction may not be 
closely enough aligned to the streamwise direction because of wing-induced spanwise velocity 
components. Convergence instabilities, especially in areas with large cell-aspect-ratios in the i-j 
computational plane and/or around the grid singularity along the wing leading edge extension 
outboard of the wing tip, may exist. To correct this situation additional time-like damping in the 
spanwise or r] direction is required. This is easily and efficiently accomplished with the following N 
operator modification 


a 


NC" jk = - a-8 i R l 


a- 5 v Sj 5 V \\a' - S(T k 8( - a z E ^ 


c n 


(33) 


where the a parameter in the third factor has been changed to a ' which is defined to be 


a' = a+a d a v {£,T},Q 

i.e., the sum of the original a parameter and a new quantity, a d a v . The quantity a v is a three- 
dimensional distribution array determined by trial and error that is permanently set, and a d is a user- 
specified constant that allows the time-like dissipation level to be increased or decreased. The 
quantity a d a v effectively provides additional time-like dissipation to the iteration process and 
stabilizes many of the difficulties encountered for swept wing computations. It also serves to stabilize 
the leading edge region of the computation in the vicinity of the local iteration and the region 
outboard of the wing tip along k= 1 where the grid cell aspect ratios are very large. The increment in 
computational cost for this modification is negligible. 


§Several versions of the leading edge local iteration algorithm have been implemented and produce 
about the same results in terms of convergence efficiency. For brevity only one is presented. 
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NUMERICAL RESULTS 

Problem setup and arid arrangement 

To evaluate the attributes of the transonic full potential algorithm just presented the familiar ONERA 
M6 wing geometry is chosen. As mentioned in the previous section, a C-H grid topology is utilized; C 
topology in the chordwise direction and H topology in the spanwise direction. All grids have been 
generated using the HYPGEN grid generation code 14 which uses a fast hyperbolic grid generation 
scheme described in Steger and Rizk 15 and Chan and Steger. 16 Best grid results are typically 
achieved (especially for the coarse grid cases) by generating finer grids than required and then taking 
only every second or third point as desired. This is especially important for the normal-like direction in 
order to keep grid lines from crossing during the grid marching process. For all of the results 
presented herein the grids were generated using triple the desired number of points in the normal- 
like direction and then reduced in size after the grid generation was completed by taking every third 
point. In addition, because of the large amount of numerical dissipation required for the grid 
generation process when marching large distances from the initial data surface, grid skewness 
develops in the grid downstream of the wing trailing edge. To partially alleviate this situation, the grid 
point distribution downstream of the wing trailing edge along each ^-coordinate line is redistributed 
using cubic spline interpolation. Despite this extra computational work, grid generation times range 
from a few seconds to a few tens of seconds on a Cray C-90 single processor computer. 

Figures 1-5 show selected views of a typical grid about the ONERA M6 wing generated using 
HYPGEN. The grid displayed is a relatively fine grid consisting of 452,925 total points, with 225 points 
in the wrap-around direction, 61 in the spanwise direction, and 33 in the normal-like direction. Figure 
1 shows a blowup of the wing surface grid in planform including a portion of the grid outboard of the 
wing tip and downstream of the wing trailing edge. In this figure only every third point in the wrap- 
around direction and every second point in the spanwise direction are displayed. Clustering in the 
wrap-around direction at the wing leading and trailing edges and in the spanwise direction at the wing 
tip are clearly visible. Figure 2 shows a close up view of the grid along a constant-chord surface near 
mid-chord. In this figure only every second point in the spanwise direction and every second point in 
the normal-like direction are displayed. This view shows the approximate treatment at the wing tip. 
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selected grid lines in the normal-like and spanwise directions (225X61X33=452,925 points). Only 
every second point in the spanwise direction and every second point in the normal-like direction are 
displayed. 

Figures 3 and 4 show two views of the wing in the wing-root cross-sectional plane; Fig. 3 shows a 
blow-up of the wing cross-section with surrounding grid and Fig. 4 shows an extreme blow-up of the 
same view. For clarity only every third point in the wrap-around direction and every second point in 
the normal-like direction are displayed in Fig. 3 while all points are displayed in Fig. 4. The grid is not 
orthogonal (especially downstream of the trailing edge) for the reason mentioned above. For all cases 
presented in this study the outer grid boundary is placed approximately 12 root chords above, below, 
upstream, and downstream of the wing and approximately 5.5 root chords outboard of the wing tip. 
An outer boundary position study conducted for the present approach and presented in Ref. 17, 
showed these distances to be acceptable for typical transonic wing computations. 

The normal and streamwise tangential grid spacing on the wing surface are approximately controlled in 
the present grid generation approach. For the grid of Fig. 4 (as well as all grids used in the present 
study) the streamwise tangential spacings at the wing leading and trailing edges are set to 0.3 and 0.5 
times the average tangential spacing, respectively. A smooth distribution of tangential grid spacing 
values in between the leading and trailing edges is obtained by numerically solving a fourth-order bi- 
harmonic differential equation for the tangential grid distribution. Use of such an approach allows easy 
specification of not only the location but also the spacing at both endpoints of this distribution. The 
resulting distribution is then interpolated onto the wing surface using cubic spline interpolation. The 
normal spacing all around the wing surface is set equal to the average tangential spacing. This value is 
automatically reduced somewhat at the wing leading edge to better match the streamwise grid 
clustering and to account for the inevitable flow gradients that exist in this area. Thus, the cells all 
around the wing surface are (on the average) approximately square. 

Figure 5 shows a similar view as that of Fig. 4 except the constant-span station is outboard of the wing 
tip. Thus, the grid singularity outboard of the wing tip at the wing-extension leading edge is clearly 
visible. Across the k= 1 boundary at this span station flow periodicity is applied, i.e., all flow variables 
on the wing extension surface at grid point / must equal the corresponding flow variables on this 
surface at NI-i+1. 
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Fig. 3 Blow-up view of the wing root grid showing grid-cell clustering around the wing surface 
(225X61X33=452,925 points). Only every third point in the wrap-around direction and every second 
point in the normal-like direction are displayed. 



Fig. 4 Extreme blow-up view of the wing root grid showing grid-cell clustering around the wing surface 
(225X61X33=452,925 points). All grid points are displayed. 
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Fig. 5 Extreme blow-up view of a typical chordwise cross-sectional grid outboard of the wing tip 
showing grid-cell clustering around the wing extension and the grid singularity at the wing-extension 
leading edge (225X61X33=452,925 points). All grid points are displayed. 


One last point about the grid displayed in Figs. 1 -5 is in order. These grid views, as well as the Mach 
number contour plots to be shown subsequently, have been generated using PLOT3D 18 . As such, 
a special Euler-like O-vector solution must be generated from the velocity potential solution. This is 
straight-forwardly accomplished from the stored values of density and from special velocity 
component values computed from the velocity potential. Because the density values are computed 
and stored at half points in the -coordinate direction (i.e., at i+1/2,j,k ), the O-vector values are also 
computed and stored at these same half points. Using half points in the £ -coordinate direction is 
also important because this allows use of the smallest computational stencil in the streamwise 
direction, which reduces smearing caused by the Q-vector post-processing computation. Likewise, 
the PLOT3D grid coordinate file is recomputed using simple second-order averages and is stored at 
these same half points in the £ direction so as to be consistent with the O-vector file. Thus, the grids 
displayed in Figs. 1-5 have been moved a small distance from where they have actually been used. 
This is why a grid space exists at the leading edge in Figs. 3-5 instead of a grid line as required by the 
leading edge flow-solver iteration scheme described in Eq. 32. 


Hybrid and Second-Order Scheme Comparisons 

In this section differences in accuracy between the hybrid spatial differencing scheme and the 
second-order-accurate spatial differencing scheme are examined. Keep in mind that the hybrid 
scheme is a second-order-accurate, centrally-differenced scheme at all subsonic points and a first- 
order-accurate, upwind-differenced scheme at all supersonic points. The first-order aspect in 
supersonic regions of flow is exclusively tied to the <£; -coordinate direction. The 77 - and £ -coordinate 
directions are always second-order accurate and centrally differenced in both subsonic and 
supersonic regions of flow. For typical weak-shock transonic applications (amenable to simulation 
using the full potential formulation), the number of supersonic points is usually around 2-4% of the 
total number of points. Thus, even though the hybrid scheme does have first-order-accurate regions, 
most of the solution behaves as if it were second-order accurate. Nevertheless, differences caused 
by these "first-order pockets of flow" can be important and are now discussed. 


Chordwise pressure coefficient distributions over the ONERA M6 Wing for both the hybrid and 
second-order-accurate schemes (as well as experimental results from Ref .19) are compared in Fig. 6. 
These comparisons are presented at four semi-span stations: 2y/b=0.2 0.44, 0.65, and 0.95. Both 
computed results utilize the grid presented in Figs. 1-5. The flow conditions for this standard case 
( m — 0.84, a = 3.06°), produce a rather benign flow field with a moderate-to-weak shock wave 

svstem with’ only a slight amount of shock/boundary-layer interaction. As can be seen in Fig. 6 the two 
numerical solutions are in excellent agreement at all subsonic locations of the flow i.e., on the entire 
lower wing surface and on the upper wing surface at and aft of the normal shock wave. In these 
reqions the computational/experimental agreement is also generally quite good, except at the normal 
shock wave where some disagreement exists primarily due to viscous effects.§ Agreement between 
the two computational results in the supersonic region (as expected) is not as good. G en e ral| y- 
forward shock, which is a supersonic-to-supersonic oblique shock, is "smeared out by the hybrid 
scheme as exemplified in Figs. 6b and 6c. The second-order scheme does a better job in this region 
capturinq the weak oblique shock wave in about the right location in Fig. 6b and about 6-7 /o of chord 
downstream of the experimental location in Fig. 6c. This second-order oblique shock location in Fig. 
6c (2y/b=0.65) is only about 2% of chord downstream of the Euler results of Refs. 20-24. 

Fffflct of C. Co. and (X o n Solution Accuracy 

In this section the effect of various parameters associated with the hybrid and second-order spatial- 
discretization schemes is examined. These parameters are essentially fixed from case to case, and 
thus are not meant to change as user-specified coefficients. However, it is of interest to see the 
effect their variation has on the solution to better appreciate the differences between the hybrid and 
second-order schemes. In all comparisons the 2 y/b = 0.65 station from the ONERA M6 Wing at 
Moo = 0.84, a = 3.06° will be used in this evaluation. This is an appropriate station because the 

sensitive supersonic-to-supersonic shock is a prominent feature at this location on the wing. 

Fiqure 7 shows a series of pressure coefficient comparisons for the hybrid scheme involving different 
values of the parameter C [defined by Eq. (10)]. Larger values of C produce larger amounts of first- 
order dissipation in the supersonic region of flow and smaller values produce smaller amounts o 
dissipation. Thus, as expected the largest value of C(C=1.7) completely smears the supersomc-to- 
supersonic shock wave and significantly "rounds out" the top of the normal shock. Conversely, the 
lowest value of C (C= 0.5) does a somewhat reasonable job at the supersonic-to-supersonic shock but 
produces a rather significant overshoot at the normal shock. In addition, stability for the lowest value 
of C is marginal requiring about 1000 iterations for convergence (more than triple the number of 
iterations required for solution convergence using the default value of C=0.9). Values of C smaller 
than 0.5 are unstable for this case. 


§Most of the computational/experimental disagreement at the normal shock wave (both in terms of 
shock position and strength) is caused by viscous effects. This qualitative conclusion is reached by 
comparing the present results with a variety of Euler results 20 ' 24 . In these Euler-based computational 
studies the normal shock position and strength are in between the full potential and experimental 
results but are more closely in agreement with the full potential results. That is, both types of invi^id 
computations produce shocks that are too strong and too far downstream relative to experiment. This 
overprediction of shock strength and position is a typical characteristic of transonic inv'scid me hods 
when used for computations that contain shock/boundary-layer viscous effects. The fact that the full 
potential results produce shocks which are somewhat stronger and further downstream of the Euler 
results implies that this additional (although generally smaller) error is a direct result of the 
isentropic/irrotational assumptions that are inherent in the full potential formulation. 
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AXIAL DISTANCE. X/C 

b) 2y/b=0.44 




Fig. 6 Pressure coefficient comparisons at selected semi-span wing stations showing differences 
between the hybrid and second-order schemes, ONERA M6 Wing, M 00 = 0.84, a = 3.06°. Both 
computations utilized the same grid consisting of 225X61X33=452,925 points. 

Figure 8 shows a series of pressure coefficient comparisons for the second-order scheme also 
involving different values of the parameter C. The supersonic domain pressure variations are not as 
dramatic for this case as compared to that of the hybrid scheme (Fig. 7). This is due to the small effect 
that the C parameter has in the second-order scheme, i.e., the term it multiplies is much smaller. 
Nevertheless, for smaller values of C the overshoot problem at the normal shock still exists. In 
addition, smaller values of C produce less stable convergence, and in particular, C= 0.5 diverged for 
this case. For second-order cases the default value of C is 1 .3. This improves convergence stability 
somewhat with a small sacrifice in accuracy. 

Figures 9 and 10 present a similar set of pressure coefficient comparisons for the second-order 
scheme involving C2 variations for Fig. 9 and C3 variations for Fig. 10. These parameters are 
inherently associated with the second-order scheme only and are defined by Eqs. (12) and (14), 
respectively. The effect of C£ on the supersonic domain pressure distribution as seen from Fig. 9 is 
quite small. Any value, including a very large value, produces acceptable results. The default for this 
parameter is 1.0. The effect of C3 on the supersonic domain pressure distribution as seen from Fig. 
10 is more pronounced. The smallest value (C3=1.0), which eliminates any growth in the second- 
order dissipative term into the flow field interior [see Eq. (14)], produces a very strong supersonic-to- 
supersonic shock and an oscillatory solution between the two shocks. Slightly elevated values of C3 
eliminate this behavior. The default value of C3 is 1 .08. 
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Fig. 7 Pressure coefficient comparisons at the 
2y/i>=0.65 semi-span wing station showing the 
effect of Con solution accuracy (hybrid scheme), 
ONERA M6 Wing, = 0.84, a = 3.06°, grid 
consists of 225X61X33=452,925 points. 



Fig. 8 Pressure coefficient comparisons at the 
2y//b=0.65 semi-span wing station showing the 
effect of Con solution accuracy (second-order 
scheme), (^=1.0, C 3 =1.08, ONERA M6 Wing, 
= 0.84, a = 3.06° , grid consists of 

225X61X33=452,925 points. 




Fig. 9 Pressure coefficient comparisons at the 
2y/i>=0.65 semi-span wing station showing the 
effect of Oi on solution accuracy (second-order 
scheme), C=1.3, C 3 =1.08, ONERA M6 Wing, 
= 0.84, a = 3.06°, grid consists of 
225X61X33=452,925 points. 


Fig. 10 Pressure coefficient comparisons at the 
2y/fc>=0.65 semi-span wing station showing the 
effect of on solution accuracy (second-order 
scheme), 0=1.3, C^=1.0, ONERA M6 Wing, 
= 0.84, a = 3.06°, grid consists of 
225X61X33=452,925 points. 


Grid Refine ment Study 

A arid refinement study for the ONERA M6 wing geometry utilized in the previous section is 
described next. Four grids, outlined in Table 1 , are used for this purpose. In Table 1 , Nl is the number 
of arid points in the wrap-around direction, NJ is the number of points in the spanwise direction, and 
NK is the number of points in the normal-like direction. The total number of points ranges from just 
under 60 000 to over 875,000. The number of points defining the wing surface ranges from 1445 to 
8651. The second-finest grid from Table 1 (case L4) is identical to the grid used in the last section 
and to the grid presented in Figs. 1-5.. 
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Table 1 . Summary of grid statistics for the grid refinement study . 


CASE 

GRID DiMENS 

ONS 

TOT SURF 

TOTAL 

NO. 

— 

NJ 

NK 

POINTS 

POINTS 

L2 

113 

31 

17 

1445 

59551 

L3 

169 

46 

25 

3175 

1 94350 

L4 

225 

61 

33 

5577 

452925 

L5 

281 

76 

41 

8651 

875596 


Selected chordwise pressure coefficient distributions are presented in Fig. 1 1 at four different semi- 
span locations, 2y/b=0.20, 0.44, 0.65 and 0.95. In each plot solutions for all four grids described in 
Table 1 are compared with experimental results from Ref. 19. As can be seen from Fig. 1 1 most of the 
solution variation affected by grid refinement is associated with the supersonic flow domain (all values 
of -C p above 0.3269). The shocks get sharper and the pressure minimums are more well defined as 
the grid is refined. The largest solution variation associated with the grid refinement process is at 
2y/b=0.95. This is probably a direct result of the approximate treatment of the wing tip geometry. The 
entire wing tip region between 2y/b=0.95 and 1.00 is represented by less than two grid cell widths in 
the span direction for the coarsest grid. Thus, as the grid is refined in this relatively high gradient 
region the solution improves remarkably. 





AXIAL DISTANCE. X/C 

c) 2y/b=0.65 



Fig. 1 1 Pressure coefficient comparisons at selected semi-span wing stations showing the effect of 
grid refinement (second-order scheme); ONERA M6 Wing, M = 0.84, a = 3.06°. 
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Another more qualitative technique for showing the effects of grid refinement on solution accuracy is 
shown in Fig. 12 where Mach number contours are displayed for the ONERA M6 Wing upper surface. 
These contours were generated using the PLOT3D plotting program 18 with contours plotted in 
0 025 increments. A contour plot is displayed for each grid described in Table. 1. Figure 12a shows 
contours for the coarsest grid, Fig. 12b the second-coarsest grid, Fig. 12c the second-finest grid, and 
Fig. 12d shows contours for the finest grid. Evolution of the solution with grid refinement, in 
particular, the upper-surface shock system, is clearly evident from this sequence of plots. The 
supersonic-to-supersonic oblique shock is the most noticeable feature that forms and sharpens as 
the grid is refined. More on grid refinement will be presented in the last results section where lift and 
drag coefficient variations with grid refinement will be compared with Euler results. 



Fig. 12 Mach number contours on the upper surface of the ONERA M6 Wing showing the effect of 
grid refinement (second-order scheme), AM = 0.025; = 0.84, a = 3.06°. 
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Convergence Efficiency 


Convergence history results are presented next showing the spatial discretization scheme's effect on 
convergence efficiency. In addition, the time-dissipation coefficient a ^ [see Eq. (33)], which affects 
convergence stability and efficiency, but not the spatial accuracy, will also be explored. Figure 1 3 
shows a plot of maximum residual, average residual and the lift coefficient as a function of iteration 
number for both the hybrid and second-order spatial differencing schemes. Both of these results 
used the L4 grid described in Table 1. In addition, the solution parameters [mainly a i , see Eq. (26)] 
have been optimized by a trail-and-error process for both of these solutions. Each symbol in Fig. 1 3 
represents 16 iterations in the convergence history for each solution. This corresponds to two 
complete applications of the a sequence, which for all cases presented herein contained eight 
elements, i.e., M=8. 



UJ 

o 


Fig. 13 Convergence history comparisons between the hybrid and second-order schemes, ONERA 
M6 Wing, = 0.84, a = 3.06°, grid consists of 225X61X33=452,925 points. 

As can be seen from Fig. 1 3, the average residual is consistently about two orders of magnitude 
below the maximum residual. This characteristic allows a simple check for healthy convergence. 
When one or more points in the solution fail to converge smoothly (perhaps due to a localized grid 
problem) the difference between the maximum and average residuals will increase. The location of 
the maximum residual at this point in the convergence will then identify the region of difficulty. The 
average residual history curve is also somewhat smoother than the maximum residual curve, which (of 
course) is due to the high-frequency error content associated with the maximum residual. In fact, if 
values were plotted for every iteration, the average residual curve would have about the same 
smoothness, but the maximum residual curve would be much more oscillatory, peaking for each value 
of at as a consequence of high-frequency error growth. This high-frequency error growth is 
accompanied by a substantial reduction in the low-frequency error, and thus, is quite important for 
overall fast convergence. The lift convergence history (also shown in Fig. 13) epitomizes how the 
convergence process eliminates the low-frequency error content in the solution. Lift changes are 
small for values of a near a h , but are quite large for values of a near a i . Generally, a typical lift 
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convergence history curve from the present algorithm involves very rapid initial growth, a moderate 
overshoot and then a highly damped oscillation to the final answer. 

Convergence histories for both solutions displayed in Fig. 13 posses the general characteristics 
described above. However, the hybrid scheme converges about 40 to 100% faster than the second- 
order scheme. This behavior is typical and is associated with increased high-frequency error 
components that exist in the second-order scheme because of reduced dissipation. For the specific 
convergence histories displayed in Fig. 13 the lift converges to within ±0.1% (a suitable condition for 
plotable accuracy) in 128 iterations for the hybrid scheme and in 240 iterations for the second-order 
scheme. These iteration counts correspond to approximately 23 and 43 sec of computer time, 
respectively, on a single processor Cray C-90. Again, for the convergence histories given in Fig. 1 3, 
the maximum residual is first reduced below a value of 10' 7 (an alternate convergence criteria that 
represents tighter convergence than plottable accuracy) in 352 iterations for the hybrid scheme and 
in 496 iterations for the second-order scheme. Thus, depending on the convergence criteria (and 
many others can be used as well), the hybrid scheme is between 1.4 to 1.87 times faster than the 
second-order scheme for the case presented in Fig. 1 3. 

The last topic for this section is to investigate and determine what effect the a d parameter has on 
solution convergence. In this regard, Fig. 14 shows several convergence histories for cases with 
different values of <x d ranging from 0.4 to 3.0. All other solution parameters are fixed at default 
values. As can be seen a d =0.4 produces an unstable result almost immediately. The other three 
values (0.8, 2.0, 3.0) all produce stable results with similar trends. Clearly, a d = 3.0 produces slower 
convergence than any other value. The two middle values of oc d (0.8 and 2.0) produce similar 
convergence histories with a d = 0.8 being somewhat faster for tighter convergence levels (especially 
if the maximum residual is used in the convergence criteria), and a d = 2.0 produces faster 
convergence for less tight levels of convergence. In terms of lift convergence, cc d = 2.0 is the fastest 
producing convergence to 0.1% of the final lift value in 368 iterations; a d =0.8 requires 384 
iterations; and a d =3.0 requires 480 iterations. 



o 100 200 300 400 500 

ITERATION NUMBER, n 


Fig. 14 Convergence history comparisons showing the effect of a d on convergence (second-order 
scheme), ONERA M6 Wing, = 0.84, a = 3.06°, grid consists of 225X61X33=452,925 points, A) 
a d = 0.4, B) a d = 0. 8 , C) a d =2.0,D) a d =3.0. 
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Lift and drag comparisons with Euler formulations 


This section presents the effects of grid refinement on lift and drag for the standard case under study, 
i.e., the ONERA M6 Wing at M «, = 0.84, a = 3.06°. In addition to results from the present full potential 
formulation, a number of other results collected from the literature are also included. The other-result 
formulations include both Euler and full potential methods using both structured and unstructured 
grids (Refs. 8, 20-25). In addition, there are newly computed full potential results (both hybrid and 
second-order) utilizing a multi-zone chimera variation of the present approach. In this approach the 
flow domain is divided into two grid zones: an inner C-H-topology grid surrounding the wing and an 
outer Cartesian-like grid that connects the inner grid to the far field using chimera interpolative 
boundary conditions. See Ref. 17 for more information about the full potential chimera grid scheme 
used in computing these results. The grids used in this series of full potential chimera computations 
are described in Table 2. The inner grid dimensions are given by Nil, NJ1, NK1, and the outer grid 
dimensions are given by NI2, NJ2, NK2. For each two-zone grid in Table 2 there is a single-zone grid 
in Table 1 with an exactly matching surface grid, both in terms of number and distribution of points. 
Thus, the surface solutions for these two sets of computations are comparable, despite the fact that 
the interior grid distribution/topology is quite different. 


Table 2. Summary of arid statistics for the two-zone chimera grid refinement study 


fgjf 

I 


MENS 


Wm 


G2L2 

101 

23 

9 

1445 

39185 


37 

19 

26 



G2L3 

151 

32 

13 

3175 

115066 


55 

25 

38 



G2L4 

201 

41 

17 

5577 

260547 


73 

33 

50 



G2L5 

251 

50 

21 

8651 

494872 


91 

41 

62 




The variation of lift and drag with grid refinement for the above list of methods are displayed in Figs. 1 5 
and 16, respectively. Values plotted along the horizontal axis are determined by taking the total 
number of wing surface nodes to the -1/2 power, which yields an approximate value for the average 
surface-cell grid spacing. Since surface quantities are being compared, it was decided to use this 
technique for computing the grid cell size. In this way structured and unstructured grid results can be 
compared? All results associated with open symbols utilize structured approaches, and all results 
associated with closed symbols utilize unstructured approaches. Results from Refs. 20-25 are Euler 
and from Ref. 8 (as well as the present results) are full potential. As can be seen from Fig. 15 the 
scatter is quite large. In fact, there is no clear cut single value for the asymptotic lift coefficient 
established by Fig. 15. A "best guess" range for the asymptotic lift coefficient limit is between 0.292 
and 0.300. In actuality, the full potential formulation, because of its inherently isentropic nature which 
leads to the prediction of slightly stronger shock waves, should produce an asymptotic lift value that is 
slightly larger than the Euler value. However, it is difficult to reach such a conclusion from Fig. 15. 


§ It is difficult to determine whether this technique for measuring grid refinement is appropriate or not 
for unstructured grid results or whether unstructured and structured grid results can be compared on 
this basis. Nevertheless, with this disclaimer, these comparisons are presented for the reader to 
evaluate. 
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AVERAGE SURFACE CELL SPACING, (N SURF )' 1/2 

Fig. 15 Variation of lift coefficient with grid refinement for a variety of different approaches including 
both Euler and full potential methods, ONERA M6 Wing, M ^ =0.84, a = 3. 06°. Open symbols are 
structured methods and closed symbols are unstructured methods. 



0 0.01 0.02 0.03 0.04 0.05 0.06 

AVERAGE SURFACE CELL SPACING, (N SURF )''' 2 

Fia. 16 Variation of inviscid drag coefficient with grid refinement for a variety of different approaches 
including both Euler and full potential methods, ONERA M6 Wing, =0.84, a = 3.06°. Open 
symbols are structured methods and closed symbols are unstructured methods. 
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The variation of drag§ with grid refinement for the above list of methods is displayed in Fig. 16. The 
same set of conventions as used in Fig. 15 are retained in Fig. 16. Again, the scatter is quite large; 
even larger than for the lift. The largest single-method variations in drag are from the unstructured grid 
results, which in some cases produce relatively large errors on moderately fine grids. For the drag 
results of Fig. 16, establishment of a definitive asymptotic value is even more difficult than for the lift. 
Nevertheless, a "best guess" range for the asymptotic drag coefficient limit is between 90 and 1 1 5 
counts. The full potential asymptotic drag value, because of its stronger shock, should be slightly 
larger than the Euler value for this case. It is interesting to note that most of the individual full potential 
drag values are below most of the Euler values, but that the trends do seem to put the full potential 
drag asymptotic limit above the corresponding Euler value as expected. The two-zone full potential 
approach (G2 curves) produces generally more favorable drag results [relative to the single-grid 
results (G1 curves)] because the inner grid generation, also performed with HYPGEN, produces a 
better quality grid near the wing surface in the chimera approach.§§ 

A more quantitative picture of the present full potential lift, drag, and pitching moment results, 
including the newly computed two-zone chimera results generated with the Ref. 17 algorithm, is 
presented in Tables 3 and 4. Also included is a variety of convergence statistics for each case. No 
attempt has been made to optimize convergence for these cases (except for the L4 grid cases 
already presented in Fig. 13). All computer timings are from a single processor of the Cray C-90 
computer using level 3 optimization for the inline, vector, and scalar options of the cf77 compiler. The 
scheme used for each computation is listed in the "SCH" column. An “H" corresponds to the hybrid 
scheme and a "2" corresponds to the second-order scheme. For the two-zone cases the first symbol 
corresponds to the scheme used for the inner grid zone and the second symbol corresponds to the 
scheme used for the outer Cartesian-like grid zone. Besides the lift, drag, and pitching moment 
information displayed in Tables 3 and 4, there are five additional quantities: the number of iterations 
required to achieve lift to within ±0.1% of the final value (nun - ), the computer time required for nuFT 

iterations (tuFT). the number of iterations required to achieve an average residual level of 10' 8 
(hravg). the computer time required for nRAVG iterations (tRAVG). and the computer time required for 
solution overhead (to/H)- The overhead time corresponds to grid generation, solution initialization 
including metric computation, and PLOT3D input file generation (a more difficult task for a potential 
solver than an Euler solver). The overhead time is larger for the single-zone results because the grid 
generation is more expensive. In the single-zone case the entire grid is generated using HYPGEN, 
and then (as mentioned previously) all points downstream of the wing trailing edge are redistributed to 
improve grid quality. In the two-zone case HYPGEN is used for only the inner grid zone (generally 
fewer than half the total points) without any point redistribution and the other Cartesian-like grid zone 
is generated using a very fast algebraic approach. 

A definite degradation in solution convergence efficiency for the second-order scheme relative to the 
hybrid scheme can be seen in Tables 3 and 4 for most levels of grid refinement. This is true no matter 
which convergence criteria is chosen as was already established for the L4 grid in Fig. 13. The two- 
zone second-order results show less degradation relative to the two-zone hybrid results. This is 
primarily due to the fact that the outer grid zone in all of these cases actually used the hybrid spatial 
scheme. Using a two-zone chimera grid arrangement with the second-order scheme for the inner grid 
and the hybrid scheme for the outer grid is attractive because it is both efficient and accurate at least at 
the wing surface where the solution is essentially second order. Convergence efficiency 
comparisons between the two-zone and the single-zone results are mixed. Sometimes the single- 
zone results converge faster and sometimes the two-zone results are faster. One observation from 
this section is the speed with which the present full potential results are obtained. Solutions times 


§ It should be stressed that the drag being discussed here is inviscid drag consisting of wave drag and 
induced drag. No drag components due to viscous effects are included. 

§§ The grid quality near the wing surface for the two-zone chimera approach is better because HYPGEN 
doesn't have nearly as far to march to the outer boundary (one chord instead of 12). Thus, even 
though there are fewer normal-direction points in the two-zone inner grid relative to the single-zone 
grid, the stretching is less which contributes to better grid quality. In addition, because the marching 
distance is much smaller, less numerical smoothing for the grid generation process is required, again 
contributing to better grid quality. 
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range from 2 sec to less than 3 min. On a typical L4 grid for a two-zone case, each solution is obtained 
in 30-40 sec of computer time. 


Table 3 Summary of computational statistics for the present single-grid full potential results displayed 
in Figs. 15 and 16, ONERA M6 Wing, AC = 0.84, a = 3.06°. Computer times are from a single 


CASE 

NO. 

SCH 

c L 

C D 

Cm 

nLIFT 

tuFT 

(sec) 

HRAVG 

tRAVG 

(sec) 

tO/H 

(sec) 

L2 

H 

.2822 

.0140 

-.1656 

64 

2 

128 

4 

7 

L3 

H 

.2877 

.0120 

-.1687 

160 

14 

192 

17 

17 

L4 

H 

.2905 

.0114 

-.1704 

128 

22 

176 

30 

39 

L5 

H 

.2922 

.0113 

-.1716 

272 

92 

256 

86 

60 

L2 

2 

.2864 

.0148 

-.1677 

112 

4 

176 

6 

7 

L3 

2 

.2911 

.0125 

-.1699 

224 

19 

256 

22 

17 

L4 

2 

.2930 

.0118 

-.1711 

240 

41 

320 

54 

39 

L5 

2 

.2943 

.0115 

-.1721 

480 

162 

432 

146 

60 


Table 4 Summary of computational statistics for the two-zone, chimera-grid full potential results 
displayed in Figs. 15 and 16, ONERA M6 Wing, = 0.84, a = 3.06°. These results use the chimera 


CASE 

NO. 

SCH 

Cl 

CD 

Cm 

nLIFT 

tLIFT 

(sec) 

ORAVG 

tRAVG 

(sec) 

tO/H 

(sec) 

G2L2 

H/H 

.2839 

.0103 

-.1653 

96 

3 

224 

8 

2 

G2L3 

H/H 

.2880 

.0106 

-.1683 

160 

11 

176 

13 

5 

G2L4 

H/H 

.2910 

.0108 

-.1705 

208 

26 

224 

28 

10 

G2L5 

H/H 

.2923 

.0109 

-.1715 

288 

66 

256 

58 

i is 

G2L2 

2/H 

.2892 

.0113 

-.1678 

128 

4 

176 

6 

2 

G2L3 

2/H 

.2916 

.0111 

-.1696 

160 

11 

224 

16 

5 

G2L4 

2/H 

.2937 

.0112 

-.1713 

272 

34 

288 

36 

10 

G2L5 

2/H 

.2944 

.0112 

-.1720 

384 

87 

400 

91 
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CONCLUDING REMARKS 

In conclusion, a new scheme for solving the full potential equation has been presented and evaluated 
using a standard three-dimensional transonic wing computation. The new scheme includes both a 
hybrid spatial discretization option which is second-order accurate in subsonic regions of flow and 
first-order accurate in supersonic regions and an option which is fully second-order accurate 
irregardless of flow type. The new second-order scheme utilizes a solution limiter, somewhat similar to 
Euler flux limiters, to maintain stable operation at shock waves and other solution extrema. 

The new iteration algorithm utilized in the present study is a variation of the AF2 scheme especially 
designed for obtaining transonic wing solutions on C-topology grids. Its special design removes (or at 
least controls) a stability limitation that exists in other implementations of this scheme. A key feature of 
the present scheme is the coupling of two different AF2-scheme variations, one above the wing and 
the other below the wing, using a local iteration at the wing leading edge. A typical transonic wing 
solution involving 450,000 point grids can be obtained in as little as 22 sec of computer time on a 
single processor of a Cray C-90 computer. 

Numerous numerical results are presented including a grid refinement study to fully demonstrate the 
new scheme's capabilities. The results indicate that the new algorithm is a viable technique for solving 
the full potential equation and could provide a very fast computational tool for the aerodynamic 
analysis and design of geometrically complex aerodynamic configurations. 
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RECOMMENDATIONS FOR FURTHER STUDY 


As stated in the introduction, the primary motivation for undertaking the present research path is to 
establish a full potential chimera capability, suitable for obtaining aerodynamic data for complex 
shapes. Thus, extending the present scheme's geometric handling capability is one area that needs 
additional work. In particular, development of a general methodology that is capable of handling 
multiple lifting surfaces with their associated vortex sheets is first on the list. 

Another area of additional work is associated with viscous effects. In order to make the present 
capability suitable for realistic aerodynamic analysis, the inviscid flow assumption must be removed. 
Thus, a boundary layer correction capability must be added to the present methodology. 

Application of the present capability to a variety of problems ranging from fast, nonlinear wind tunnel 
wall corrections to complex-geometry design applications is the next logical step for this research. 

Finally, improvements in computational efficiency are required. Although not discussed in this report, 
the processing rate of the present algorithm is only marginally adequate, being about 300 MFLOPS 
on a single processor of a Cray C-90 computer when a fine grid is used. A variety of efforts in the 
parallel processing area could lead to dramatic improvement in this processing rate and a 
correspondingly dramatic improvement in the utility of this technology. 
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